arXiv:1505.00073vl [cs.GR] 1 May 2015 


BIJECTIVE DEFORMATIONS IN R N VIA INTEGRAL CURVE COORDINATES 


1 


Bijective Deformations in M" via Integral Curve 
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Abstract —We introduce Integral Curve Coordinates, which identify each point in a bounded domain with a parameter along an integral 
curve of the gradient of a function / on that domain; suitable functions have exactly one critical point, a maximum, in the domain, and 
the gradient of the function on the boundary points inward. Because every integral curve intersects the boundary exactly once, Integral 
Curve Coordinates provide a natural bijective mapping from one domain to another given a bijection of the boundary. Our approach can 
be applied to shapes in any dimension, provided that the boundary of the shape (or cage) is topologically equivalent to an n-sphere. 
We present a simple algorithm for generating a suitable function space for / in any dimension. We demonstrate our approach in 2D 
and describe a practical (simple and robust) algorithm for tracing integral curves on a (piecewise-linear) triangulated regular grid. 

Index Terms —Deformation, bijection, homeomorphism, coordinates, integral curves 
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1 Introduction 

Shape deformation is a widely studied problem in com¬ 
puter graphics, with applications in animation, finite 
element simulation, parameterization, interactive model¬ 
ing, and image editing. Deformations are useful for 2D 
and 3D shapes, as well as higher-dimensional shapes. 
Deforming the animation of a solid model is in fact a 4D 
problem, as the animating 3D solid is itself a 4D shape. 
Interpolating parameterized modeling spaces, such as 
the space of human body shapes 111, can be an arbitrary 
dimensional problem. 

In the version of the problem that we study, a bound¬ 
ary or "cage" is created around a shape. As this cage 
is manipulated, the interior is also deformed. The cage 
may be identical to the shape's boundary, which has one 
fewer dimension than the shape itself, and is typically 
more convenient, as the cage may be simpler (fewer 
vertices) or be free of undesirable properties (such as 
a non-manifold mesh or high topological genus). 

One important yet elusive property for deformations 
is bijectivity. A bijective function is a one-to-one map¬ 
ping such that every point in the undeformed figure is 
mapped to a point in the deformed figure and vice-versa. 
Bijectivity ensures that the shape never flips "inside-out" 
as a result of deformation. 

We propose a technique that creates a guaranteed 
bijective deformation of the shape given a bijective de¬ 
formation of its boundary or cage, in any dimension. To 
do so, we introduce Integral Curve Coordinates, which 
identify each point in the shape by its position along an 
integral curve of the gradient of a suitable function /. 
Suitable functions are Lipschitz continuous, have only 
a single critical point, a maximum, in the interior, and 
have gradients that point inward along the boundary of 
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the domain (or cage). For such functions, integral curves 
never meet unless they are identical, and every integral 
curve can be uniquely identified with a point on the 
boundary. Given a bijective deformation of the boundary 
points, we construct a suitable function / and trace the 
integral curves through the interior points to extend the 
bijective deformation into the interior of the shape. Our 
approach can be viewed as a generalization of Xia et al. 
|2j from star-shaped domains to arbitrary contractible 
domains. 

Our contributions are: 

• Integral Curve Coordinates which naturally pro¬ 
duce a guaranteed bijective deformation between 
two contractible domains in M n , given a bijective 
deformation of their boundaries. 

• A simple algorithm to create a scalar function space 
on a regularly discretized contractible domain in 
M n ; functions in this space have exactly one critical 
point, a maximum. 

• A practical (simple and robust) algorithm for tracing 
2D integral curves on a (piecewise linear) triangu¬ 
lated regular grid. 

We note that the bijections produced by Integral Curve 
Coordinates, while correct, are not fair. The recent work 
of Schuller et al. |3| and Aigerman and Lipman |4| 
complements our own by improving the fairness of a 
deformation in a manner that preserves bijectivity. 

2 Related Work 

Deformation is a well-studied problem in computer 
graphics. One approach to shape deformation warps 
the entire ambient space, which induces a deformation 
for the coordinates of the shape @-0 In contrast, 
intrinsic shape deformation approaches operate in terms 
of relative coordinates, without regard for the ambi¬ 
ent space j8|-[10|. Cage-based deformations enclose the 
shape in a sort of structural boundary (which may 
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simply be the shape's actual boundary). This boundary, 
when deformed, induces a deformation of the interior 
points m —1151. Several approaches allow a combination 
of cages, line segments, and sparse points while still 
computing an intrinsic shape deformation 1161, [17]. 

One approach to cage or boundary deformation in¬ 
volves the generalization of barycentric coordinates 
[18 [-(26). Traditionally defined, barycentric coordinates 
express points in the interior of a simplex (triangle in 
2D, tetrahedron in 3D) as the weighted average of the 
vertices of the simplex. The simplex can be thought of 
as a simple cage around its interior; the deformations 
induced by modifying vertices of the simplex are bijec- 
tive (so long as the simplex remains non-degenerate). 
Generalized barycentric coordinates extend this idea to 
more complex polygons or polyhedra than triangles and 
tetrahedra, but cannot guarantee bijectivity (27). 

2.1 Bijectivity 

Few deformation approaches guarantee bijectivity. Bijec¬ 
tivity guarantees that the deformation does not invert 
(flip inside-out) or collapse any part of a shape. In 
two-dimensions, Weber et al. [ [24] introduced conformal 
mappings based on complex coordinates. Xu et al. [ |28] 
introduced a thorough solution to the challenging prob¬ 
lem of determining vertex locations for the interior of 
a triangular graph, so that no triangles flip inside-out. 
This is akin to solving for a bijective deformation with 
the additional constraint that mesh connectivity remains 
unchanged. 

While generalized barycentric coordinates cannot 
guarantee bijectivity (27), Schneider et al. [29] recently 
introduced a technique that splits a deformation into a 
finite number of steps, each of which is implemented 
via generalized barycentric coordinates. In practice, the 
technique creates bijective mappings at pixel accuracy. 
One limitation, however, is that a continuous non-self- 
intersecting interpolation between the source and target 
cages is required, which is an unsolved problem for 
dimensions greater than two. A theoretical analysis of 
this technique has not yet been performed in 3D (or 
higher). 

Lipman's restricted functional space, introduced ini¬ 
tially in 2D [30] and extended to 3D and higher by 
Aigerman and Lipman |4) can be used to find bijective 
deformations of piecewise linear meshes similar to a 
desired one. The technique projects an input simplicial 
map onto the space of bounded-distortion simplicial 
maps. While they do not prove that their iterative al¬ 
gorithm convergences in general, if it does converge, it 
guarantees a bijective deformation. 

The approach of Schuller et al. [3] prevents inverted 
elements (i.e. guarantees local injectivity) by augmenting 
any variational deformation approach with a non-linear 
penalty term that goes to infinity as elements collapse. 
The approach therefore requires a continuous deforma¬ 
tion from an initial shape, and is incompatible with hard 
constraints (such as a required target pose). 
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Several approaches have been presented for creating 
bijective texture maps with positional constraints [31]— 

m 

The approaches of Weinkauf et al. [341 and Jacobson 
et al. [35] both create restricted function spaces that 
prevent undesirable extrema (maxima and minima), but 
cannot control the placement of saddles. Jacobson et 
al. applied their technique to shape deformation. Our 
approach creates a restricted function space with exactly 
one extrema and no saddles, in any dimension. 

The goal of our work is to introduce a cage-based 
deformation technique that guarantees bijectivity in any 
dimension, and provide a practical piecewise linear im¬ 
plementation in 2D. 

3 Background 

Definition 1: A function h: X —Y is called bijective if 
and only if it is one-to-one: 

h(xi) = h(x 2 ) xi = x 2 

and onto: 

\/y E Y, 3x e X such that h(x) = y. 

Bijective functions always have an inverse. In the 
literature, the desired bijective deformations are actually 
homeomorphisms, which are bijective functions h such 
that both h and hr 1 are continuous. 

For deformations, the one-to-one property (injectivity) 
is the most difficult to achieve. A lack of injectivity 
manifests as regions of X that collapse or invert (flip 
"inside-out") when mapped via h. For piecewise linear 
shapes, this corresponds to collapsed or inverted ele¬ 
ments (triangles in 2D, tetrahedra in 3D, etc). Expressed 
symbolically, a function h must have a Jacobian whose 
determinant is everywhere positive to be injective: 

det{J h (x )) >0 ViEl (1) 

Where the determinant is negative, h has locally flipped 
X inside out (inverted). Where the determinant is zero, 
h has locally collapsed X , though not inverted it. (In 2D, 
a triangle collapses to a line or point.) 

For piecewise linear shapes, the Jacobian is piecewise 
constant within each element and can be constructed as 
follows. Consider a /c-dimensional simplex P (triangle 
in 2D or tetrahedron in 3D). P has k + 1 vertices 
vq,vi , Ignoring the overall translation of h, the 

mapping h(P) can be expressed via matrix multipli¬ 
cation. Choose a vertex of P arbitrarily (vo, without 
loss of generality). The matrix we seek maps — vo 
to h(vi) — h(v o). Assuming that P is non-degenerate, 
vi — vo ,..., Vk — vo are linearly independent. The matrix 
M such that M(vi — ^o) = h{vf) — h(v o) is, therefore, the 
product of two matrices whose columns are vertices: 

M = [h(v i) - h(v 0 )\h(v 2 ) - h(v 0 )| • • • | h(v k ) - h(v 0 )] 
x [t>i - v 0 \v 2 - «o| • • • \vk ~ vo] -1 
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The determinant of M determines whether h is locally 
injective within P. If the determinant is zero, h has 
locally collapsed P to, in 2D, a line segment or point. 
If the determinant is negative, h has inverted P. 

4 Overview 





Function Creation 


Integral Curve Tracing Output 


Our goal is to deform the interior of a shape given a 
deformation of its boundary or enclosing cage. Formally, 
our approach takes as input 

1) The boundary dD of a contractible domain D in 

M n . 

2) A homeomorphism h deforming dD. 
and extends h to the interior of D. 

To accomplish this, we introduce Integral Curve Co¬ 
ordinates (formally defined below) that identify a point 
p in the domain D with an integral curve x(t) and a 
parameter t p such that x(t p ) = p. The integral curves 
follow the gradient of a special function f D : ^x(t) = 
V/£>(x(t)). The function f D is constructed such that 
has exactly one critical point in D , a maximum, and 
its gradient on the boundary points inward. Integral 
curves of the gradient of a function are sometimes called 
integral lines p6) . The integral curves of a Lipschitz 
continuous function are well-defined everywhere except 
critical points, and two integral curves never meet unless 
they are identical. Because every integral curve of V/d 
traces a path from a point on the boundary of the domain 
to the maximum, we can uniquely identify every integral 
curve with the boundary point it passes through. We 
denote the integral curve passing through boundary 
point b m as Xb m (t). The Integral Curve Coordinate of 
a point p is 


Tfo(p) 


(b m,t) ifp^po 

0 if p = Po 


(3) 


where x bm (£) = p and p 0 is the maximum point of fn- 

Since distinct integral curves never meet, X is a bi- 
jection from points in the domain to Integral Curve 
Coordinates. Given a bijective deformation h of the 
boundary of D , h(dD) = dD' , we can similarly construct 
a special function f D i on D' to create 2/ d , . Transforming 
from Xf D to Xf D , can therefore be achieved by mapping 
the to /i(6 m ). 

To summarize, the steps to extend h to a point p in 
the interior of D are as follows (Figure [lj: 

1) Create a function f D : D M with a single critical 
point, a maximum, and whose gradients on the 
boundary point inward. 

2) Create a similar function f' D for D'. 

3) Compute the Integral Curve Coordinate I p = 
Ifni P) by tracing the integral curve of V fo in both 
directions from p. 

4) Transform I p to V, by mapping the boundary point 
b m of the Integral Curve Coordinate with h (or, if 
p = po, identifying the unique maximum po of 
with the unique maximum p' 0 of f' D ). 


Fig. 1. An overview of our approach. Given a shape D 
and a homeomorphism of its boundary ft: dD dD', we 
first create two functions, one inside dD and one inside 
dD', each with only a single critical point, a maximum. 
To find the new position of a point p e D, we trace 
the integral curve passing through p up to the maximum 
and down to a boundary point b m . We then trace the 
integral curve from h(b m ) on the boundary of D' to the 
maximum. The new position of p is the point located the 
same fraction along the integral curve in D'. Repeating 
this process everywhere allows us to extend the boundary 
homeomorphism h to the interior of dD'. 


5) Compute Xj * (I ' p ,) by tracing the integral curve of 
V/d' from the boundary point to the maximum. 
(To warp points backwards from D' to points in D , swap 
D with D' and h with ft -1 in the above steps.) 

Our implementation of this approach, described in the 
following sections, creates functions on a regular grid 
discretization of the domain. We trace integral curves on 
a piecewise linear interpolation of the function values. In 
the piecewise linear setting, gradients are piecewise con¬ 
stant and discontinuous across piece boundaries. This 
guarantees monotonic interpolation (no spurious critical 
points within an element) and greatly reduces numerical 
issues in tracing integral lines. However, it violates the 
assumption that integral lines never meet, and can result 
in regions of the shape collapsing (though not inverting). 
In M 2 , any of a number of C 1 -continuous monotonic 
interpolation techniques would eliminate this problem 
(37)-(40). 

5 Function Creation 

The first step in our deformation approach is the cre¬ 
ation of a suitable function upon which to trace inte¬ 
gral curves. Namely, we require functions with a single 
critical point, a maximum, and whose gradients on the 
boundary point inward. Our discrete implementation 
creates a function space satisfying the above require¬ 
ments. The function space takes the form of a set of in¬ 
equality relationships on edges of a regular grid enclos¬ 
ing the shape boundary or cage. We note that harmonic 
functions, which obey the strong maximum principle, 
are suitable in 2D, but not in higher dimensions as they 
may contain spurious saddle points. 

Our approach first chooses a grid vertex to be the 
maximum (Section \5.1) and then generates inequality 
relationships and vertex values such that all other points 
are regular (Section [5^2| . We then trace integral lines on 
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a piecewise linear interpolation of the function values 
(Section [6). 

5.1 Maximum Selection 

Any grid vertex can be selected to be the maximum. 
However, because integral lines converge towards a 
maximum, we wish to choose a maximum vertex in 
the center of the shape. We find such a point with 
the grassfire transform [411 To do so, we construct a 
regular grid covering the mesh. The grassfire transform 
iteratively "burns away" the boundary of a shape; we 
use it to find the farthest point from the boundary. Any 
vertex in the final iteration can be chosen arbitrarily. 

5.2 Assigning Function Values 

Once the maximum is chosen, we assign function values 
to all grid vertices inside and just enclosing the shape 
boundary or cage. To do so, we build a function space 
containing functions that have a single critical point, a 
maximum, and gradients on the boundary that point 
inward. The function space takes the form of a cousin 
tree of inequality relationships on edges of the regular 
grid. 

Definition 2: Let T be a tree defined on a subset of a 
regular grid. We call T a cousin tree if every pair of axis- 
aligned neighbor vertices in the tree are related to each 
other as parent-child or as tree cousins (Figure |2|. Two 
vertices are tree cousins if they are neighbors in the grid 
and their tree parents are also neighbors in the grid. 

tree parent 

© O — O 


© ©— © 
tree child tree cousins 

Fig. 2. Cousin tree vertex relationships: Two vertices u,v 
that are axis-aligned neighbors in the grid must be either 
parent-child (left) or tree cousins (right). Tree cousins’ 
parent vertices must also be axis-aligned neighbors in the 
grid. 

Our algorithm for creating a cousin tree, illustrated 
in Figure]3j is a form of breadth-first search (BFS) that, 
for all vertices on the frontier in a given iteration, always 
expands first along the first coordinate axis, second along 
the second coordinate axis, and so on. (The order of 
coordinate axes is not important so long as it is consistent 
across all iterations.) Pseudocode is as follows: 

1: procedure CreateCousinTree (maximum Vert ex) 

2: > Store the tree as a set of directed edges: 

3: edges <— {} 

4: frontier <— {maximumVertex} 

5: functionValue <— empty dictionary 

6: functionValue[maximumVertex] <— 1 
7: while frontier is not empty do 
8 : frontierNext <— {} 


9: for dimension in fixed dimension order do 

10: for v in frontier do 

11: candidates <— { unvisited neighbors of v along 

dimension direction } 

12: edges <— edges U {(v,n) \ n G candidates} 

13: filtered <— {n E candidates if n within boundary} 

14: frontierNext <— frontierNext U filtered 

15: > Optional: directly assign function values: 

16: function Value [ filtered] function Value [v] — 1 

17: end for 

18: end for 

19: frontier <— frontierNext 

20: end while 

21: return edges, functionValue 

22: end procedure 

When a vertex expands, it becomes the tree parent 
of the vertices it expands into. In 2D, this means that 
all potential children along the x-axis are connected to 
the tree in one round. The next round, all potential 
children along the y -axis are connected to the tree, if 
they have not already been included. While vertices 
outside the boundary are not added to the frontier of 
the breadth-first search, they are added to the tree and 
considered children of all adjacent vertices within the 
domain; this ensures that vertices outside the boundary 
have smaller values than vertices inside, which cre¬ 
ates inward-pointing gradients at the boundary. See the 
supplemental materials for a proof that this algorithm 
indeed creates a cousin tree. 






Fig. 3. A cousin tree is created via breadth-first search 
expanding outwards from a maximum (red circle). 


The cousin tree represents our function space; the edge 
from a parent to a child vertex represents the inequality 
f (parent) > f (child). Note that the Manhattan distance 
lies within this function space, as breadth-first search 
can be used to directly compute the graph distance to 
a given node. See the appendix for a proof that this 
function space contains no critical points other than the 
maximum. 

One the cousin tree is created, we can assign values 
to grid vertices such that the tree parent always has 
greater value than the tree child. We have experimented 
with various approaches to assigning function values: 
directly assigning the L\ (tree) distance, and solutions 
to the Laplace and bi-Laplace equations subject to vari¬ 
ous boundary conditions and the cousin tree inequality 
constraints. These approaches are evaluated in Section [7] 


6 Tracing Integral Lines 

With a suitable function fo in hand, we are ready to trace 
integral lines to convert a point p in the interior of the 
shape to its Integral Curve Coordinate Xf D ( p) = (b m ,t), 
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where b m is the point on the boundary reached by 
tracing the integral line through p in the decreasing 
direction (gradient descent), and t is the fraction of arc- 
length along the integral line from b m to the maximum, 
reached by gradient ascent from either p or b m . (Note 
that given a piecewise linear boundary or cage, we 
store b m in barycentric coordinates with respect to the 
boundary piece it belongs to.) 

Since our function values are provided on a regular 
grid, we require monotonic interpolation of the values 
within grid cells, so that spurious critical points are 
not introduced by the interpolation. For our 2D imple¬ 
mentation, we perform piecewise linear interpolation via 
a regular triangulation that divides each grid cell into 
two triangles. Notably, piecewise linear interpolation is 
monotonic and has constant gradient, which mitigates 
problems resulting from numerical accuracy when trac¬ 
ing integral lines and simplifies arc-length parameteri¬ 
zation. 

Because the gradient is constant within each triangle, 
one could trace the integral line from a given starting 
point by simply computing the intersection of the ray 
from the starting point in the gradient direction with the 
boundary of the triangle. However, naive implementa¬ 
tion of this numerical approach may create discrepancies 
and inconsistencies such as loops due to limitations in 
floating point precision. In order to obviate numerical 
precision issues, we separate the topological calculation 
(which edge of the triangle the integral line passes 
through) from the geometric calculation (the coordinates 
of the point on the edge). Our topological calculations 
are based on comparing values and computing the signs 
of cross products, which are amenable to symbolic per¬ 
turbation schemes j42}-{44) for robust, exact evaluation)^] 

Our algorithm traces an integral line by iteratively 
computing the sequence of triangle edges it intersects 
(and points on those edges). As a special, initial case, 
the integral line may originate from a point inside a 
triangle, but thereafter will be tracked via the triangle 
edges it intersects. Our algorithm proceeds as follows, 
and is illustrated in Figure |4] Given a point p\ located 
on the edge e Bs of a triangle B, we determine the 
next edge by computing the sign of the cross product 
between (a) the gradient of the function within B and 
(b) the vector from p 1 towards the vertex of the triangle 
opposite the edge e Bs (the dotted red line in Figure [3] 
right). The sign of the cross product of these two vectors 
determines the outgoing edge of the integral line (Table [l] 
with edge labeling given by Figure [5|. This robustly and 
stably computes the outgoing edge. The point on the 
outgoing edge (p 2 in Figure [4] right) is then determined 
via standard line /line intersection computation (which 
is made simple since triangles edges are aligned with 
grid edges or diagonal). Tracing proceeds in the triangle 
on the other side of the outgoing edge. 


1. We did not implement such a scheme, as we did not encounter 
numerical issues with our floating point implementation. 


In the special, initial case of an integral line originating 
at a point p 0 inside a triangle A, the sign of the cross 
product of the gradient is computed with each of the 
three vectors from the point to the triangle's vertices 
(Figure |4j left). The sign of the cross products determines 
which vectors the gradient is to the left and right of, 
which indicates the outgoing edge accordingly. The point 
on the edge can then be computed numerically, and the 
general case of the algorithm proceeds. 


&Aj 


eAi 




eBi 


Fig. 4. An integral line is traced (left) from p 0 inside 
triangle A. The gradient direction (the purple arrow) is 
compared to the red dotted lines, indicating which edge 
of A is intersected. Line/line intersection tells us the 
numerical location of the next point p 1 lying on a triangle 
edge, which is the general case. Tracing then proceeds 
inside the opposite triangle B (right), and the gradient 
direction only needs to be compared to one dotted red 
line. 



Fig. 5. The edge labeling for triangles in the regular grid, 
such that the lookup table in Table [j] can be used to 
compute the outgoing edge when tracing integral curves. 

Finally, because our piecewise linear domain is not 
C 1 , integral lines may meet. This manifests in our inte¬ 
gral line tracing algorithm as a triangle whose gradient 
points backwards against the incoming edge. We call 
such edges compression edges. Without special care, the 
integral line would zig-zag or staircase between the two 
triangles. We detect this by testing whether the triangle's 
gradient points towards the incoming edge (e.g. right 
back out of the triangle). If so, then tracing return to the 
face opposite the incoming edge, and the next point on 
the integral line is the endpoint of the edge pointed to¬ 
wards by the gradient. Because our triangles triangulate 
a regular grid, the dot product of the gradient with the 
edge itself is a simple sign test or comparison between 
components of the gradient. This procedure simulates 
the integral curve bouncing between the two triangles as 
their gradients merge, skipping the intermediate steps. 
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Edge Sign Next Edge 


1 + 3 

1 - 2 

2 + 1 

2 - 3 

3 + 2 

3 - 1 


TABLE 1 

Our algorithm for tracing integral lines jumps from one 
triangle edge to the next (Figurej4| right). With edges 
labeled according to Figure[5pthe next edge is 
determined by the sign of the cross product between the 
function gradient and the vector towards the triangle’s 
opposite corner. 


7 Evaluation 

We experimented with various functions for Integral 
Curve Coordinates. Figure [6] depicts the computed func¬ 
tion, integral lines, and deformations for the following 
functions subject to cousin tree constraints: L\ (tree) 
distance and solutions to the Laplace and bi-Laplace 
equation with boundaries equal to zero; and solutions 
to the Laplace equation with boundaries equal to zero, 
without constraints imposed by the cousin tree. We wish 
for integral curves to be maximally separated. Solutions 
to the Laplace equation when boundary values are 
fixed to zero result in integral curves that intersect the 
boundary orthogonally, resulting in the most separation 
between integral lines. There is little difference between 
solutions to the Laplace equation with and without the 
cousin tree constraints. 

Figures [7jj9] compare the results of our deformation 
approach to Complex Barycentric Coordinates |24|, Con¬ 
trollable Conformal Maps 1141, Composite Mean Value 
Mappings (29}, and Locally Injective Mappings (3). Our 
deformation, while less fair, can be generalized to any 
dimension. Although some integral lines do flatten out 
as a result of deformation, no inverted elements are 
created by our approach. Our deformation can be used 
to compute a correct starting configuration for recent 
techniques which improve the fairness of deformations 
while preserving bijectivity 

The long and narrow shape in Figure [9] emphasizes 
the degree to which the location of the maximum affects 
the overall deformation. We believe that replacing the 
single maximum point by a skeleton or medial axis is a 
fruitful direction for future research. 

Performance Our experiments were written in unop¬ 
timized Python and executed on a 2 Ghz Intel Core 
i7. In all examples, we compute function values on a 
50-by-50 discrete grid. Performance is dominated by 
per-pixel integral curves tracing. Our examples contain, 
on average, 38,850 pixels, and took approximately 30 
minutes each. There are large performance improve¬ 
ments to be obtained by implementing integral curve 
tracing in a compiled language and by parallelization. 
Further performance improvements could be obtained 


bi-Laplace Laplace Laplace 

equation equation equation 

with with without 

Li (tree) cousin tree cousin tree cousin tree 
undeformed distance constraints constraints constraints 


a* &* $ + $ + $ + 






Fig. 6. Integral Curve Coordinates computed using var¬ 
ious functions. For each object, the first row depicts the 
uv deformation map in the red and green channels; the 
second row warps a checkerboard pattern; and the third 
row displays the function values (background lightness), 
maximum location (red circle) integral lines (blue), the 
cousin tree, and expansion edges (purple, Section [H}. 




by deforming all points along an integral curve at once, 
rather than wastefully re-tracing integral curves for each 
point along it. Finally, an approach based on advecting 
the boundary could efficiently trace all integral curves 
at once. 

8 Limitations 

In our piecewise linear implementation, compression 
edges (Section [6} occur when the gradients on either side 
of an edge point towards the edge. Similarly, we call an 
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Fig. 7. An extreme deformation of a square, a counter¬ 
example to the bijectivity of generalized barycentric sys¬ 
tems (45j. Left column: A square shape with a checker¬ 
board pattern, Integral Curve Coordinates with a Laplace 
equation not subject to cousin tree constraints. Middle col¬ 
umn: Harmonic coordinates [46] , Controllable Conformal 
Maps |[14J. Right column: Locally Injective Mapping using 
Dirichlet and Laplacian energy (3). Only Integral Curve 
Coordinates and Controllable Conformal Maps generate 
a bijective deformation. 



Fig. 8. Left: The undeformed shape. Center: A de¬ 
formation computed using Complex Barycentric Coordi¬ 
nates (24). Right: A deformation computed using Integral 
Curve Coordinates; the function is the solution to the 
Laplace equation with the boundary values constrained 
to zero, subject to the cousin tree constraints. Complex 
Barycentric Coordinates produce very smooth deforma¬ 
tions, but are limited to 2D. 


edge an expansion edge when the gradients on either side 
point away from it. When tracing integral lines downhill, 
multiple integral lines converge at an expansion edge. 
This leads to multiple integral lines intersecting the same 
boundary point; the points along these integral lines 
will "collapse" as a result of the deformation. Expansion 
edges rarely occur for the L\ function values or the 
Laplace equation, except near small concavities on the 
boundary. Expansion edges are visualized in Figure [6] as 
purple edges, and typically only occur near concavities 
of the boundary. With C 1 or G 1 function interpolation, 
expansion edges (and compression edges) would no 
longer occur. A looser requirement than continuity is 
simply that the gradients on either side of an edge never 
point away from each other; a tangible solution for this 
relaxed condition is unclear. 

Our grid discretization of the boundary or cage may 


is 

Fig. 9. From left to right: The undeformed shape, a defor¬ 
mation computed using Integral Curve Coordinates with 
a Laplace equation not subject to cousin tree constraints, 
and the deformation computed using Composite Mean 
Value Mapping [29]|. 



lead to problematic boundary gradients near sharp an¬ 
gles (< 45°). One solution is to warp space with a 
simple "plaid deformation" such that (a) grid vertices 
are positioned exactly at boundary vertices and (b) sharp 
angles are non-uniformly scaled and eliminated. 

9 Conclusions and Future Work 

Integral Curve Coordinates provide a new approach 
for bijective shape deformation based on tracing the 
integral curves of functions with one critical point, a 
maximum. While the deformations produced by our 
approach are not as fair as, for example. Controllable 
Conformal Maps [14), they are bijective in all dimen¬ 
sions. The fairness of our deformations can be improved 
in a bijectivity-preserving manner via the recent, com¬ 
plementary work of Schuller et al. [3| and Aigerman 
and Lipman (5). We believe that fairer functions may 
be found in our function space. One approach may be 
to compute a compatible skeleton for the deformed and 
undeformed shapes and treat the entire skeleton as the 
maximum. 

Our approach is restricted to cage-based or boundary 
deformations. In the future, we would like to extend our 
approach to other control structures, such as points and 
bone skeletons, which are intuitive to manipulate and 
can have far fewer vertices than a cage. One could trace 
integral curves of a smoothed distance functions from 
the control geometry (47). 

Our piecewise linear implementation, while prevent¬ 
ing inverted elements ( det(M ) < 0 in Equation 2), does 
not prevent collapsed elements ( det(M ) = 0), which are 
caused by compression edges. We would like to address 
this in the future with C 1 or G 1 monotonic interpolation 
of functions values, or by simulating the infinitessimal 
separation of integral curves (48) and then perturbing 
the resulting deformation to correct collapsed elements. 

We would also like to explore modifications to the 
cousin tree constraints. The constraints we compute, 
while correct, are not unique. Thus, we envisage an 
iterative procedure that updates the constraints and the 
function values. Jacobson (27) explored such an iterative 
constraint modification scheme in an analogous setting 
to good effect. A similar iterative scheme may also 
remove expansion edges. 
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Finally, we would like to implement our approach in 
higher dimensions, applying it to problems such as the 
animation of volumetric models. 
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Appendix 

In the following, the domain is a regular n-dimensional 
grid with edges parallel to the X 2 , £ 3 , • • •, x n axes. We 
call a subdomain ball-like if the volume enclosed by its 
cells is contractible. 

Theorem 1: Let D be a subdomain that is ball-like. 
Let T be an associated cousin tree. Assume all nodes 
on the boundary of the subdomain have value greater 
than all grid neighbors outside the ball (including along 
imaginary diagonal grid edges). Then by rooting T at 
the maximum and assigning monotonically decreasing 
(unique) values to internal nodes of T along the path 
from the root (preserving the values at the root and 
the leaf), there can be no grid maximum, minimum, 
or saddles at any tree node with 4- or 8 -connectivity. 
(In dimensions > 2, we assume monotonic interpolation 
within a hypercube; i.e. 2 n-connectivity.) 

Lemma 1: A node v of cousin tree T on subdomain D 
with parent in the +Xi grid direction has a child in the 
—Xi grid direction, unless the neighbor in the — Xi grid 
direction is outside D. 

Proof: Without loss of generality, assume node v has 
its tree parent u in the Ax^ grid direction. Let w be the 
node in the —x\ grid direction. Assume v is not w's tree 
parent. 

9 — 9 —? 

o—©—o 
o —® o 


The cousin rule tells us that v's tree parent u and w's 
tree parent must be grid neighbors. Yet the only grid 
neighbor of w within 1 grid edge of v's parent of u is 
v. Thus we have reached a contradiction, and v must be 
w's parent. □ 

Lemma 2: For any non-boundary tree node v, its grid 
neighbors (even with the addition of diagonal edges on 
the hypercube face planes) with greater value all belong 
to the same connected component (partitioned according 
to greater/less than v's value) and its grid neighbors 
(even with the addition of diagonal edges) with lesser 
value all belong to a second connected component. 

Proof: Without loss of generality, assume v has par¬ 
ent u in the +£2 direction. We now consider two cases, 
v is on the boundary of D and v is not on the boundary 
of D. 

Case v is not on the boundary of D: 


It follows from Lemma [l] that v has a child w in 
the — x 2 direction. As paths are monotonic, Value(u ) > 
Value[v ) > Value(w). Since we wish to prove that v 
will have exactly two connected components partitioned 
according to greater/less than v's value, we can restrict 
our examination to grid neighbors of u,v,w in the +x\ 
direction, again without loss of generality; the greater 
value connected component will have to include u and 
the lesser value connected component will have to in¬ 
clude w. We now consider two subcases, v is not the 
tree parent of its x\ -direction grid neighbor b and v is the 
tree parent of b. To clarify, when determining connected 
components, neighbors along diagonal edges may be 
considered; however, tree edges are always grid edges, 
and neighbors in any sense except when computing 
connected components are only along grid edges. 



Because of the cousin rule, b's parent must be a, since 
a is the only grid neighbor of b within 1 grid edge of 
v's parent of v. Lemma [l] tells us that c must be the tree 
child of b. Due to the monotonicity of values along tree 
paths, Value(u ) > Value(v ) > Value{w) and Value(a) > 
Valuefb ) > Value(c). The following diagram depicts all 
possible greater/less than relationships between a, 6 , c 
and v. 


©- 

—© 

©— 

—® 

®— 

- 4 ) 


As we can see, in all possibilities, there are exactly 
two connected components partitioned according to 
greater/less than v's value. The component with values 
greater obviously contains v's tree parent u, and the 
component with values lesser obviously contains v's tree 
child w. 

Subcase: v is the tree parent of b. 


©— 

—® 

©— 

—® 

®— 

—© 


The cousin tree imposes the following restrictions on 
a and c. a cannot be the tree child of b as u and a would 
not be tree cousins. Since b cannot be a's tree parent or 
its tree child, b and a must be tree cousins. Therefore 
a's tree parent must be a grid neighbor of v . The only 
possibility is u. Since c cannot be the tree parent of b 
(resp. w), it must be the tree child or cousin of b (resp. 
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w). Therefore c must be the tree child of one and the 
tree cousin of the other. In either case, the monotonicity 
of values along tree paths implies Value(v) > Value(c) 
as well as Value(v ) > Value(b ), Value{v) > Value{w) f 
and Value(u) > Value(v). The value of a may be greater 
than or less than v, but in both cases there are exactly 
two connected components (partitioned according to 
greater/less than v's value). As in the earlier subcase, 
the component with values greater obviously contains 
v's tree parent u, and the component with values lesser 
obviously contains v's tree child w. 

Case v is on the boundary of D: 

If v's grid neighbor w in the — x 2 direction is in D, 
then Lemma [l] applies and v must be the tree parent of 
w and hence Value(v) > Value(w). Otherwise, w is not 
in D, but by assumption Value{v ) > Value(w). Thus v's 
neighbor in the +x 2 direction has value greater than v, 
and v's neighbor in the — x 2 direction has value less than 
v. Following the same argument as in Case v is not on the 
boundary of D, we can again restrict our examination to 
grid neighbors of u,v,w in the +x 2 direction (without 
loss of generality). 



—® 


—® 


A) 


There are 33 valid cousin tree configurations among 
the 2 4 possible boundary conditions given a,b,c,w can 
each be outside D. They are presented in Figure |To| 

In every case there are exactly two connected com¬ 
ponents partitioned according to greater/less than v's 
value. The component with values greater obviously 
contains v's tree parent u, and the component with 
values lesser obviously w. □ 

Proof: Lemma [2] tells us that any internal tree node 
v has one connected component of nodes with value 
greater than v's and another connected component with 
values lesser (0 "folds"). It follows directly then that v 
cannot be a minimum, maximum, or saddle. □ 

Theorem 2: A breadth-first search (BFS) in a ball-like 
domain D where at each BFS generation all x\, then all 
x 2/ then all x 3 , and so on, grid edges are explored in 
order constructs a cousin tree. 

Proof: By induction. Let n represent the number of 
BFS generations of growth from the root node. After step 
n, nodes reached at BFS generation n — 1 have all of 
their grid neighbors in the BFS tree or outside D; grid 
neighbors in the BFS tree will be shown to satisfy the 
cousin tree definition. 

After step n = 1, the root node has become the parent 
of all grid neighbors. The root node satisfies the cousin 
tree constraints by being the BFS parent of all grid 
neighbors. 

Assume true for n = k. Nodes reached at BFS gener¬ 
ation k — 1 have all grid neighbors also in the BFS tree 


with cousin tree definition satisfied or outside D for BFS 
generation < k — 1 nodes. 

We wish to show that after step n = k + 1 all nodes 
reached at BFS generation k now also have all their 
grid neighbors in the BFS tree satisfying the cousin tree 
definition or outside D. 

Consider a grid node v reached at BFS generation k. 
Let u be the BFS parent of v . u was thus reached at BFS 
generation k — 1. Consider the following arbitrary axis- 
aligned figure of v: 


o e> o 

M© - *£) S 

O 0 O 


After step n = k + 1 nodes a,c,w are also in the 
BFS tree. We aim to show that a, c, w are either the BFS 
children of v, tree cousins of v in the BFS tree, or outside 
D. 

If w (or a, c) is outside D, then it is of no concern to us. 
If w (or a, c) is inside D, then it must have been reached 
at step k — 1 , k, or k + 1 . w (or a,c) cannot have been 
reached before step k — 1 since it would have been able 
to reach v at step k — 1. w (or a, c) must be reached before 
step fc + 2 since v can reach it at step k - hi. 

Suppose node w was reached at generation k — 1. 

k „©—►£)—^ 

The inductive hypothesis tells us that w has cousin 
tree relationships to all its grid neighbors, including v. 
Yet w is not the BFS parent, BFS child, or BFS cousin of 
v (since the BFS parent of w cannot be a grid neighbor 
of u). This contradicts w having been reached at BFS 
generation k — 1. 

Suppose node w was reached at generation k. 

,0 -*£) 0 

This too is impossible since v,w are neighbors and 
cannot have the same taxicab distance to the BFS root in 
our ball-like domain D. 

This leaves us with w reached at generation k + 1. 



We first consider p as the BFS parent of w. Then p must 
have been reached at BFS generation k. Furthermore, if 
p,v were both reached at BFS generation k and w was 
reached at BFS generation fc +1, then a must be in D and 
have generation k — 1, since D is ball-like and distance 
is taxicab. But p as the BFS parent of w contradicts the 
BFS xi,X 2 ,xs ,... growth ordering implied by u as the 
BFS parent of v, since u is the BFS parent of v instead 
of a. The same argument prevents r as the BFS parent 
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Of w. Consider q as the BFS parent of w. Then q and 
v were reached at the same BFS generation while w 
between them was reached at a later BFS generation. 
This is impossible in our ball-like domain D since BFS 
generations correspond to minimal taxicab distance. The 
only remaining possibility is v as the BFS parent of w. 
This is a valid cousin tree relationship. 

Suppose node a was reached at generation k — 1. 


k-1 


© 

© dS> © 

©- *£l) ® 


The inductive hypothesis tells us that a has cousin tree 
relationships to all its grid neighbors, including v. 
Suppose node a was reached at generation k. 


k-1 



This too is impossible since v,a are neighbors and 
cannot have the same taxicab distance to the BFS root 
in our ball-like domain D. 

Suppose node a was reached at generation k + 1. 

© 


© kip © 

M©-® 

If node p is the BFS parent of a, then p was reached at 
BFS generation k. We then have the following diagram. 


®—,0—® 

„.,© —*0 -® 

But nodes v and p cannot have the same taxicab 
distance (as evidenced by their BFS generations) to the 
BFS root in ball-like D if nodes a and u differ in taxicab 
distance to the BFS root by 2. The same argument 
prevents e from being the BFS parent of a. The only 
remaining possibilities are v as the BFS parent of a and b 
as the BFS parent of a, both of which have valid cousin 
tree relationships with v. 

The analysis of the relationship between node v and 
node c follows exactly the same argument as between 
nodes v and a. 
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Fig. 10. Cousin tree configurations for a node on the 
boundary: a outside (2 cases); a, b outside; a,c outside; 
a, w outside; a,b,c outside; a,b,w outside (violates as¬ 
sumption); b outside; b, c outside (4 cases); b, c, w outside 
(4 cases); c outside (5 cases); c,w outside (5 cases); w 
outside (6 cases); a, b, c, w outside. 




